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Modern  military  technology  has  become  highly  dependent  on  a  detailed  knowledg 
of  atom-molecule  and  molecule-molecule  interactions.  This  tyoe  of  information  is 
required  in  diverse  defense  anplications  including  chemical  laser  development,  in 
the  detection  and  modeling  of  plumes,  and  in  the  decomposition  of  energetic  mater¬ 
ials.  The  description  of  forces  governing  molecular  reactions  is  provided  by 
potential  energy  surfaces.  These  surfaces  are  the  crucial  first  step  in  dynamics 
calculations  that  nrovlde  required  information  about  state-to-state  cross-sections 
and  rate  constants.  Since  potential  energy  surfaces  are  not  generally  available 
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'  ABSTRACT 

from  experiment,  the  most  reliable  approach  to  their  determination  lies. in  the 
development  and  application  of  predictive  ab  initio  quantum  mechanical  methods. 
The  following  annual  report  describes  our  research  on  the  development  of  many- 
body  perturbation  theory  (MBPT)  and  related  infinite-order  coupled-cluster  (CC) 
methods  for  notential  energy  surfaces.^ 

In  the  past  year,  several  thing^Tiave  been  accomplished.  The  couoled- 
cluster  model  has  been  generalized  to  include  all  contributions  from  single  and 
double  excitation  operators,  (CCSD),  which  means  the  model  accurately  includes 
quadruple  excitations  and  some  triple  excitation  effects.  The  other  fourth-order 
triple  excitation  terms  have  also  been  added  [CCSD+T(4)].  This  very  highly 
correlated  wavefunction  shows  excellent  promise  for  describing  energy  surfaces 
accurately.  Comparisons  made  with  full  Cl,  which  includes  all  possible  excitatior 
demonstrate  that  CCSD+T(4)  obtains  all  of  the  available  correlation  energy  surface 
for  BeH2  even  shows  that  CCSD  can  correctly  describe  highly  degenerate  bond 
breaking  processes.  ; 

In  other  work,  we  have  (1)  investigated  localized  orbital  descriDtions  of 
molecules  in  the  CCSD  model  to  begin  to  develop  approaches  that  are  more  suited 
to  large  molecules;  (2)  began  the  development  of  multi -reference  MBPT/CC  method 
to  better  describe  energy  surfaces;  and  (3)  initiated  a  study  of  infra-red  in-  j 
tensities. 

-  The  importance  of  the  infra-red  signature  to  plume  identification  and 
detection  makes  it  important  to  develop  techniques  for  accurately  computing  the 
infra-red  intensities  of  molecules,  particularly  those  of  the  transient  free- 
radical  s  that  occur  in  nlumes.  This  year  we  have  developed  and  annlied  our 
MBPT/CC  techniques  to  intensities  for  a  series  of  molecules  includinq  diatomic 
inter-halogens  and  halogen  hydrides  that  arise  in  plumes  as  components  of  spent 
fuel . 
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on 

MOLECULAR  INTERACTIONS  WITH  MANY-BODY  METHODS 

to 

AIR  FORCE  OFFICE  OF  SCIENTIFIC  RESEARCH 

from 

Quantum  Theory  Project 
University  of  Florida 
Gainesville,  Florida  32611 

I .  DESCRIPTION  OF  PROBLEM 

In  a  wide  variety  of  Air  Force  applications,  highly  detailed  information 
about  atoms,  molecules,  and  their  interactions  is  required[l-3] .  This  Infor¬ 
mation  1 s  necessary  in  problems  ranging  from  chemical  laser  development,  to  the 
detection  and  identification  of  rocket  plumes,  to  metal  clustering  and  aerosol 
formations,  and  even  to  nuclear  weapons  ef fects[l-3]. 

Probably  the  most  crucial  conponent  needed  to  understand  molecular  reactions 
is  the  potential  ener^  surfaces  that  serve  to  describe  the  attractions  among  the 
atoms  and  molecules[l J .  However,  such  information  is  not  easy  to  obtain.  A 
certain  amount  of  information  about  the  molecular  forces  near  equilibrium  in  a 
bound  molecule  is  available  from  spectroscopy.  Some  information  about  the  poten¬ 
tial  energy  surface  even  in  the  absence  of  binding  can  be  provided  from  crossed 
molecular-beam  experiments.  But,  in  general,  potential  energy  surfaces  are  not 
amendable  to  experimental  determination.  Instead,  other  types  of  experimental 
observations  such  as  kinetics  experiments,  coupled  with  very  simple  theoretical 


models  for  a  surface,  are  used  Co  Infer  pieces  of  information  abouC  Che  parame¬ 
ters  of  Che  model  such  as  whaC  Che  activation  barrier  might  be. 

The  most  direct  approach  to  obtaining  detailed  information  about  a  potential 
energy  surface  is  offered  by  predictive,  ab  initio  quantum  mechanical  calcula¬ 
tions.  However,  to  make  it  feasible  to  calculate  accurate  energy  surfaces  for 
molecules,  much  better  and  more  computationally  efficient  methods  must  still  be 
developed. 

One  such  approach,  namely  many-body  perturbation  theory  (MBPT)  [4-15]  and 
its  infinite-order  extensions  termed  coupled-cluster  (CC)  methods  [11,16-22] 
offer  a  number  of  attractive  features  that  the  more  traditional  configuration 
interaction  approaches  lack[21].  Under  AFOSR  support,  we  have  generalized  the 
theory  and  developed  very  efficient  computer  codes  to  perform  MBPT/CCM  calcula¬ 
tions  and  employed  these  for  the  first  time  in  large-scale  ab  initio  calculations 
of  potential  energy  surfaces[ 1 1,2 1] .  The  successes  of  this  effort  have  been 
substantial.  These  include  the  determination  of  a  complete  force-field  for  the 
H20  molecule,  including  all  force-constants  through  fourth-order,  that  is  suf¬ 
ficiently  accurate  that  once  improved  experiments  were  carried  out  after  our 
calculations,  many  of  the  previously  accepted  values  for  the  force  constants  were 
revised  to  be  more  consistent  with  our  predictions[22] .  Also,  a  study  of  the 

+  NH3,  and 

H^BCO  +  BH^  +  CO  was  made  that  predict  these  binding  energies  to  within  1 
kcal/mole  of  the  accepted  experimental  values  for  diborane  and  borane  car¬ 
bonyl,  and  make  a  prediction  in  the  case  of  borazane  in  the  absence  of  an 
experlment[14].  Earlier  experiments  which  gave  much  higher  values  for  the 
binding  energies  of  diborane  and  borane  carbonyl  than  we  computed  are  now 
completely  discounted.  Similar  successes  with  studies  of  the  isomerization 
energy  and  activation  barrier  of  HNC  -►  HCN[23],  and  CH^NC  -*■  CH3CN[24],  the 


binding  energies  of  the  molecules  B2H6  +  2BH3,  HjBN^  +  BHj 
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photodissociation  of  f ormaldehyde[25],  and  various  studies  of  methanol, 
methoxy,  and  the  fornyl  radical[26]  attest  to  the  reliability  of  our  MBPT/CC 
methods. 

Building  upon  this  work  supported  by  the  AFOSR,  we  carried  out  extensive 

studies  of  the  potential  energy  surface  for  the  two  Inelastic  collisions, 

3  3 

0(  P)  +  H2O  and  0(  P)  +  CO2,  under  contract  to  the  Air  Force  Rocket  Propulsion 
Laboratory,  for  the  purpose  of  obtaining  vibrational  excitation  cross-sections 
that  are  needed  in  actual  detection  devices[27]. 

Under  our  new  AFOSR  grant  starting  November  1,  1981,  we  have  initiated 
three  exciting  new  research  categories  that  represent  major  additions  to  our 
earlier  pioneering  work  with  many-body  methods  for  energy  surfaces.  These  are: 

•  Localized  orbital  descriptions  of  molecules  and  their  energy 

[' 

surfaces.  This  work  is  directed  toward  developing  the  theory 
required  for  accurate  predictions  about  large  molecules. 

•  The  development  of  the  theory  and  application  of  multi-reference 
MBPT/CC  techniques.  This  effort  is  geared  to  better  describing 
energy  surfaces  when  several  configurations  are  quite  important 
in  an  accurate  wavefunction  and  for  accurate  predictions  of 
excitation  energies. 

I 

I 

•  Predictions  of  properties  like  dipole  moments  and  polarizabilities 
and  their  derivatives  with  nuclear  displacements.  Such  derivative 
quantities  are  responsible  for  infra-red  and  Raman  intensities  and 
thereby,  offer  crucial  information  for  plume  signature  detection. 

These  properties  are  also  very  difficult  to  predict  accurately, 
placing  extreme  demands  on  the  accuracy  of  any  quantum  mechanical 
method. 


This  past  year  we  have  made  much  progress  in  each  of  these  areas.  A 
survey  of  some  specific  acconplishments  is  presented  in  Section  II,  while 
Section  III  elaborates  on  these  topics  as  well  as  describing  a  few  other  new 
results.  Section  IV  lists  the  presentations  and  publications  supported  by 
AFOSR  from  1981  and  1982.  Two  of  the  papers  we  have  published  this  year  are 
included  as  Appendices  to  describe  in  more  technical  terms  some  of  our 
research  results. 

II.  Review  of  Research  Accomplishments 

The  overall  objectives  of  this  research  program  include  the  following: 

1)  Develop  new,  more  accurate  and  more  efficient  ab  initio  quantum 

mechanical  methods  based  upon  MBPT  and  CCM  for  determining  molecular 
properties,  and  particularly,  potential  energy  surfaces  for  molecular 
interactions. 


2)  Implement  these  methods  in  highly  efficient,  transportable  computer 
codes,  to  enable  computations  on  potential  ener©  surfaces  to  be 
made  on  an  almost  routine  basis. 


3)  Apply  these  techniques  to  a  variety  of  problems  that  are  of 
interest  to  AFOSR,  and  that  serve  to  establish  the  range  of 
accuracy  for  MBPT  and  CCM  methods. 

In  line  with  these  overall  objectives,  a  number  of  accomplishments  have 
been  made  under  our  new  AFOSR  grant.  The  main  focus  for  our  effort  this 
last  year  has  been  the  development  of  formal  methods  and  new  conputer 
programs  designed  to  generalize  MBPT  and  CCM  theory.  In  particular,  the 
inclusion  of  monoexcited  clusters  in  a  full  CCSD  model  have  been  developed 
and  efficiently  implemented  within  the  CCM  mode.  Also  fourth-order  triple 
excitation  diagram  terms  including  the  diagrams  required  for  non-SCF  orbi¬ 
tals  have  been  Implemented.  Using  this  highly  accurate,  correlated  model, 
comparisons  have  been  made  with  full  Cl  (i.e.  the  best  possible  solution  in 
a  basis  set)  and  other  very  accurate  Cl  calculations  that  include  higher 


categories  of  excitations.  Since  CCSD  is  invariant  to  wide  varieties  of 
transformations,  localized  orbital  techniques  for  molecules  are  being 
investigated  to  develop  methods  for  large  molecules.  A  summary  of  achieve¬ 
ments  follows. 

A.  For  the  first  time,  the  equations  for  the  full  CCSD  model,  built 

upon  the  wavefunction,  ■  expOT^^^)  (  4>o>,  have  been  obtained 

and  Implemented  on  the  computer.  Initial  applications  were  to  J^O 
and  Bel^  where  full  Cl  solutions  were  available.  In  both  cases, 
CCSD  recovers  better  than  98%  of  the  correlation  energy  in  the 
basis  set[28],  while  the  addition  of  triple  excitations  makes  the 
agreement  almost  exact. 

B.  An  approximate  reaction  path  for  the  insertion  of  Be  into  was 
also  determined  at  the  CCSD  level  and  compared  with  full  Cl.  This 
problem  is  very  difficult  to  describe  with  a  single  reference  func¬ 
tion  since  both  Be  and  H2  have  at  least  two  important  con¬ 
figurations  along  this  path.  In  spite  of  this  difficulty, CCSD 
agrees  with  the  full  Cl  result  to  within  1  kcal/raole  (i.e.  .04 
e.v.)  for  the  full  reaction  path.  This  result  attests  to  the  sta¬ 
bility  of  the  CCSD  model  even  in  a  difficult  multi-reference 
environment! 2 9].  (See  Appendix  A  for  the  paper.) 

C.  In  the  treatment  of  large  molecules,  localized  orbitals  would 
appear  to  be  necessary  to  eliminate  the  highly  delocalized  nature 
of  the  usual  SCF  orbitals.  This  reduces  the  number  of  integrals  in 
a  calculation  down  to  a  manageable  size.  Furthermore,  such 
orbitals  are  more  conceptual  and  offer  better  prospects  for  trans¬ 
ferability  to  other  molecules.  This  year,  we  have  made  the  first 
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many- body  calculations  with  localized  orbitals,  using  the  CCSD 
model  and  easily  generated  hybrid  bond  orbitals.  Comparisons  with  SCF 
orbital  calculations  on  CH^  and  H^O  are  excellent,  with  the  loca¬ 
lized  orbital  CCSD  calculation  differing  from  the  SCF  orbital  CCSD 
result  by  less  than  1  kcal/mole[30]. 

D.  In  a  substantial  achievement,  Schaefer,  et.  al.  ,  have  carried  out  a 
Cl  calculation  involving  more  than  one  million  single,  double, 
triple  and  quadruple  excitations  on  CjH^  [31].  We  have  recently 
coapleted  a  CCSD  calculation  on  the  same  system  which  is  an  impro¬ 
vement  over  the  Cl  result  by  3  kcal/mole.  Of  this,  triple  excita¬ 
tions  account  for  about  5  kcal/mole  of  correlation  energy.  In 
addition,  CCSD  requires  1-2  orders  of  magnitude  less  computing  time 
than  the  CI[32].  (See  Appendix  B  for  the  manuscript.) 

E.  Two  full  Cl  calculations  have  recently  been  accomplished  by 
Harrison  and  Handy  [33],  one  for  H20  at  geometries  displaced  from 
equilibrium  and  one  for  the  BH  molecule.  For  these  benchmark 
results  we  have  performed  full  fourth-order  perturbation  theory 
SDTQ-MBPT(4)  and  CCSD+T(4)  where  the  fourth-order  triple  excitar- 
tions  are  added.  Our  CCSD+T(4)  results  agree  to  within  0.27 
kcal/mole  for  H20  at  equilibrium,  and  1.0  kcal/mole  for  a  displace¬ 
ment  of  one-half  the  OH  bond  length.  For  BH  in  a  larger  basis  set, 
the  difference  is  only  0.73  kcal/mole. 

F.  Preliminary  correlated  calculations  on  the  derivatives  of  the 
dipole  moment  with  respect  to  displacements  of  the  nuclei  in  H20, 
demonstrates  that  SDQ-MBPT(4)  (the  fourth-order  approximation  to 
CCSD)  can  describe  these  properties  to  within  a  mean  error  of  7%, 


III.  Summary  of  Research  Problems 

In  the  past  year  we  have  Initiated  a  localized-orbital  many-body  method 
directed  toward  large  molecules  and  we  hope,  eventually,  to  address  the 
transferability  of  parts  of  molecules  into  other  bonding  environments.  Our 
inital  efforts  have  used  very  simple  bonding  orbitals  composed  of  hybrid 
orbitals  on  the  various  atomic  centers.  In  this  way,  we  have  studied  CH^, 
CjH^,  C2H2,  C2H6>  an<*  H2°*  w^icJ>  two  papers  are  published  [30,32]  and  a 
third  one  under  preparation.  Each  of  these  molecules  have  a  well-defined 
hybridization  scheme,  which  is  why  we  chose  these  systems  for  an  initial 
study.  However,  we  have  implemented  the  Del'Re  procedure  [35]  to  define 
hybrid  orbitals  in  an  unbiased  way  regardless  of  bonding  environment  for 
more  complicated  cases. 

Our  procedure  for  making  simple  localized  bond  orbitals,  which  is  simi¬ 
lar  to  that  used  in  the  semi-empirical  PCILO  scheme  [36],  is  to  define  an 
atomic  orbital  basis  set,  hybridize  the  atomic  orbitals  to  establish  direc¬ 
tionality,  and  then  combine  the  hybrids  together  with  other  hybrid  or  ato¬ 
mic  orbitals  to  give  normal  chemical  bonds.  For  CH^,  these  bond  orbitals 
have  the  general  form 

b  =  C(sp^)  +  X  H(ls) 

where  the  weight  factor,  X  ,  reflects  the  relative  charge  distribution. 

By  doubly  occupying  the  C(ls)  core  orbitals,  and  the  four  C-H  bond  orbi¬ 
tals,  a  crude  zeroth-order  localized  reference  function  is  constructed. 
Similarly  an  antibonding  orbital  is  formed  from  the  negative  combination 


b  -  C(sp3)  -  X  H(ls) 


The  remaining  polarization  functions  In  the  basis  set  are  Incorporated  by 
suitable  orthogonallzatlon  techniques  discussed  elsewhere  [30]. 

In  the  CCSO  model,  the  wavefunctlon  Is 

^CCSD  “  exP(Ti+r2^  I  *o> 

where 

T1  “  I  t®  afl 
1 

t2  -  i  ta1yibtj 

i<J 

a<b 

Since  T^  and  T2  comnute, 

exP(T1-w;2)  |  $q>  -  exp(T2)exp(T1)  |  *Q> 

ar.d  any  single  determinant  function  can  be  Inter-related  to  any  other  <£0 
by  exp(T1),  via  Thouless'  relation  [37] 

|  lQ>  -  exp(T1)  |  $0> 

Consequently,  even  with  a  poor  choice  of  4>0  as  given  (e.g.)  by  relatively 
crude  localized  bond-orbitals,  once  the  amplitudes  (t^)  are  converged, 
one  recovers  a  good  single  reference  function,  |Q  . 

However,  If  it  Is  desirable  to  use  low-order  perturbation  methods 
instead  of  converging  to  Infinite  order  as  in  the  CCSD  model,  better  per¬ 
turbation  approximations  are  possible  when  a  better  reference  function, 

$0  ,  is  chosen.  One  way  to  do  this  while  maintaining  the  other  advantages 
of  localized  orbitals  is  to  first  do  an  SCF  calculation  and  then  localize 


the  orbitals  using  a  technique  like  that  of  Foster  and  Boys  [38]  or  Edmlston 


—  10  — 


and  Ruedenberg  [39].  This  has  Che  advantage  that  Che  computational 
reference  energy  is  optimum  since  the  localizati  .  procedure  just  amounts 
to  a  unitary  transformation  of  the  canonical  SCF  result.  Another  advantage 
is  that  this  type  of  procedure  is  completely  independent  of  geometrical 
considerations  which  make  it  easy  to  follow  a  molecule  along  a  reaction 
path  while  redefining  the  localized  orbitals  on  the  path.  We  have  not  yet 
pursued  this  localized  orbital  variation  but  we  plan  to  this  coming  year. 

Another  research  area  which  we  are  actively  pursuing  is  the  development 
of  the  multi-reference  MBPT/CC  model.  As  discussed  in  Section  II  and  Appendix 
A,  we  have  studied  the  insertion  of  Be  into  Hj  because  this  problem  has 
an  extremely  difficult  multi-reference  character.  The  Be  atom  has  a  near 
degeneracy  between  its  2s  and  2p  orbitals.  Similarly,  when  is  separated 
to  permit  the  Be  to  insert,  the  two  Hj  configurations,  lo  and  la  become 
nearly  degenerate.  Consequently,  it  is  extremely  difficult  to  see  how  one 
can  describe  this  insertion  without  including  all  of  the  important  con¬ 
figurations  from  the  separated  products  into  the  reference  sp«-ce  for  the 
problem.  For  the  molecule,  the  multi-reference  character  is  manifested  by 
the  fact  that  two  BeJ^  configurations,  la1^2a^^3a^^  and  la1^2a^23b2^ 
are  dominate  at  different  points  along  the  reaction  path  for  the  insertion, 
causing  an  ambiguity  in  which  to  choose  as  a  reference  for  any  single 
reference  approach. 

By  changing  the  reference  function  along  the  path,  CCSD  was  shown  to 
accurately  reproduce  the  full  Cl  result  for  this  problem[29].  However,  having 
to  change  reference  functions  is  undesirable,  so  for  many  problems  of  this 
type,  a  multi-reference  technique  is  requisite. 

We  have  addressed  the  multi-reference  problem  this  ye$r  in  two  ongoing 
efforts:  (1)  a  finite-order  small  reference  MBPT  approach;  and,  (2)  an 


infinite-order  large  reference  linearized  multi-configuration  coupled- 
cluster  method  (L-MCCCM).  Each  approach  has  its  advantages  and 
limitations.  Since  almost  nothing  is  known  about  the  numerical  charac¬ 
teristics  of  such  methods,  it  is  first  necessary  to  explore  alternatives. 
Besides  the  Bef^  application; we  are  also  studying  the  potential  energy  curve 
for  the  N2  molecule,  which  is  an  extremely  difficult  case.  Each  of  these 
multi-reference  approaches  also  should  enable  the  calculation  of  excitation 
energies,  as  well,  which  is  a  goal  of  this  research  program. 

In  our  work  on  properties  other  than  the  energy,  we  have  focused  on 
infra-red  intensities  due  to  their  pertinence  in  plume  detection  and  iden¬ 
tification.  Infra-red  Intensities  are  obtained  from  the  derivative  of  a  mole¬ 
cular  dipole  moment  with  respect  to  bond  and  angle  displacements  and  offer  an 
exceedingly  difficult  property  to  compute  accurately  with  even  the  best 
quantum  mechanical  methods.  Our  initial  applications  have  been  to  the 
prototype  molecule  HjO,  where,  even  for  such  a  simple  system,  the  discre¬ 
pancies  between  experiment  and  theory  have  been  severe.  In  our  calcula¬ 
tions,  SDQ-MBPT(4)  has  been  able  to  predict  these  dipole  derivatives  to 
within  a  mean  error  of  7%  while  other  calculations  have  been  in  error  by 
over  30%  [34].  We  are  currently  ivestigating  the  effect  of  infinite-order 
single  excitations  via  the  CCSD  model  on  those  predictions  and  the  effnct 
of  triple  excitations.  Other  applications  have  been  made  to  inter-halogen 
molecules  and  the  hydrogen  halides  because  of  their  relevance  to  plume 
signatures.  This  work  Includes  a  prediction  of  the  experimentally  unknown 
dipole  moment  of  the  IF  molecule. 

In  related  applications,  the  highly  unusual  HIF  inorganic  free-radical 
molecuie[ 40]  which  is  stable  to  decomposition  to  HI+F  or  IF4-H  has  been  studied 
to  locate  its  transition  state  for  the  exoergic  decomposition  to  HF+I . 
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This  very  unusual  system  has  two  nearly  equal  minima  at  an  HIF  angle  of 
127°  and  80°  suggesting  the  possibility  of  the  H  atom  tunnelling  through 
the  double-well  potential.  We  expect  to  publish  this  work  soon. 

Also,  we  have  studied  several  electronic  states  of  the  HCF  molecule, 
which  1 8  a  common  chemiluminescent  constituent  of  plumes.  A  number  of 
interesting  questions  pertain  to  the  relative  stability  of  the  singlet  and 
triplet  form  of  this  free-radical  and  the  location  of  other  low-lying 
excited  states[41].  Our  calculations  are  the  most  accurate  available  for 


this  unusual  species 
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Abstract 

The  coupled-cluster  single  and  double  excitation  model  (CCSD)  is  applied  to  an  energy  path  for  the  in¬ 
sertion  of  Be  into  Hi  and  compared  to  the  full  configuration  interaction  (FCl)  and  full  valence-multicon¬ 
figuration  self-consistent  field  (fv-mcscf)  results.  This  model  problem  is  a  severe  test  of  a  single-refer- 
encc-function  correlated  method  since  two  configurations  are  heavily  weighted  in  the  FCI  description.  CCSD 
is  demonstrated  to  describe  the  FCI  results  using  a  single  reference  function  which,  however,  changes  orbital 
characteristics  along  the  sampling  path.  In  this  case  CCSD  gives  excellent  agreement  with  the  FCI  re¬ 
sults. 


1.  Introduction 

A  number  of  applications  of  coupled-cluster  methods  (COM)  [  1  -4]  and  the  closely 
related  many-body  perturbation  theory  (mbpt)  [5-8]  for  molecules  and  their  energy 
surfaces  have  been  made  in  the  last  few  years.f"In  essentially  every  case  [10],  the 
calculation  used  a  single-reference  configuration,  usually  of  SCF  form  as  a  restricted 
(RHF)  or  unrestricted  Hartree-Fock  (UHF)  function.  With  a  single  reference  function, 
mbpt/CC.M  techniques  offer  uniform,  size-extensive,  and  highly  efficient  procedures 
for  the  treatment  of  electron  correlation. 

Although  a  number  of  approaches  to  the  multireference  mbpt/ccm  problem  exist 
[11-16],  the  applications  of  these  techniques  have  been  severely  limited  [17-19],  due 
to  the  comparative  complexity  of  the  multidimensional  approach.  As  a  consequence, 

*  Supported  by  the  U.S.  AFOSR  under  Grant  No.  82-0026. 

rFor  a  recent  review  of  CCM/MBPT  and  their  applications,  see  Ref.  9. 
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the  predominant  criticism  of  mbpt/CCM  has  been  the  practical  limitation  to  a  single 
reference  function,  particularly  in  a  region  of  bond  breaking  when  more  than  one 
configuration  is  highly  weighted  in  a  configuration-interaction  (ci)  wave  function. 
No  doubt  such  a  criticism  is  justified  in  many  cases,  however,  the  problem  deserves 
closer  investigation  to  accurately  assess  the  limitations  of  the  single-reference  approach. 
For  example,  the  coupled-cluster  doubles  (CCD)  model  has  been  shown  to  correctly 
•model  some  severely  degenerate  problems  [20]  as  has  the  coupled-cluster  single  and 
double  excitation  model  (CCSD)  [21], 

Furthermore,  the  multireference  procedure  is  not  a  universal  solution.  Besides  the 
comparative  complexity,  the  principal  objection  to  multireference  techniques  is  the 
ambiguity  in  the  choice  of  reference  configuration.  Typically,  single  a.ia  double  ex¬ 
citations  from  the  full  reference  space  will  be  chose'  to  form  the  Ct  wave  function 
However,  to  separate  a  polyatomic  molecule  like 

C 

D—A—B 

into  its  appropriate  fragments,  requires  that  the  multireference  space  contain  all 
configurations  that  are  needed  to  describe  the  component  subspecies.  Consequently, 
if  we  are  breaking  the  A — B  bond,  we  typically  need  to  include  the  configurations  that 
describe  those  fragments,  often  even  including  the  united  atom  configuration(s)  [22]. 
However,  when  we  choose  to  break  another  bond,  like  A — C  or  A — D,  we  need  a  new 
set  of  reference  configurations  which  yield  a  new  set  of  single  and  double  excitations 
to  be  included  in  the  CI.  Furthermore,  since  a  typical  objective  is  to  compare  reaction 
paths  for  the  different  channels,  the  different  levels  of  approximation  partially  in¬ 
validate  such  numerical  comparisons.  The  alternative  of  including  all  possible  reference 
configurations  for  the  plausible  reaction  paths,  which  would  eliminate  the  ambiguity, 
is  usually  not  feasible  for  polyatomic  systems  since  the  size  of  the  reference  space 
becomes  prohibitive. 

Faced  with  this  dilemma,  single-reference  approaches  which  allow  different  reaction 
pathways  to  be  treated  within  the  same  well-defined  model  should  be  examined.  After 
all,  regardless  of  the  reference  function  chosen,  the  corresponding  full  Cl  calculation 
will  give  the  ultimate  answer  within  the  basis  set.  If  one  can  introduce  the  important 
correlation  effects  by  a  suitable  treatment  of  higher  excitations,  one  should  be  able 
to  treat  even  some  very  difficult  problems  with  single-reference-function  coupled- 
cluster  techniques.  (When  degeneracies  arc  important,  the  infinite-order  ccm  can 
be  expected  to  be  more  stable  than  finite-order  mbpt  [3,8],  so  we  will  mainly  focus 
on  the  use  of  CCSD  in  this  paper.) 

We  have  previously  considered  the  BeH2  system  with  CCSD  at  a  fixed  geometry 
[21].  In  this  paper,  we  have  expanded  this  study  to  include  a  quasireaction  path  for 
the  perpendicular  insertion  of  Be  into  Hj  to  form  BeHj,  where  several  configu¬ 
rations  should  be  important  along  the  energy  surface.  This  system  is  small  enough 
that  the  full  Cl  calculation  [23,24]  can  be  performed  to  provide  rigorous  comparisons. 
In  addition,  we  report  full  valencc-multiconfiguration  self-consistent  field  (fv- 
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MCSCF)  calculations,  which  use  all  of  the  expected  important  configurations  while 
simultaneously  optimizing  the  molecular  orbitals  to  obtain  the  best  possible  variational 
fit  for  that  set  of  configurations  on  the  energy  surface  [25,26]. 


2.  Synopsis  of  CCSD 

The  CCSD  theory  is  presented  in  detail  elsewhere  [21  L^Tjiere  we  present  only  a  brief 
summary.  The  CCSD  wave  function  is 

^ccsd  =  exp(X  (  +  X2)|<^o>. 

where 

X,  = 

iM 


aU  by. 


We  make  the  convention  that  the  indices  i,j,  k,  l  represent  orbitals  occupied  in  the 
reference  configuration  0o,  while  the  indices  a,  b,  c ,  d  are  for  excited  orbitals,  with 
the  corresponding  creation  a f  and  annihilation  i  operators.  The  two-particle  amplitudes 
are  antisymmetrized 

*ob  — .  a ob  __  __  .bo 

lij  —  lji  1  ij  1  ji  » 

where  r?  are  the  one-particle  amplitudes. 

By  considering  the  projections  of  the  Schrodinger  equation  for  the  trial  wave 
function  i/'ccsd  onto  the  space  of  single  and  double  excitations,  a  set  of  nonlinear 
equations  sufficient  to  determine  amplitudes  tf  and  are  obtained.  Once  the 
equations  are  solved  for  the  amplitudes,  the  energy  is  given  by 

E  =  <0o|//exp(X,  +  X2|0o>, 


The  quantity  (ab\\ij)  =  (ai\bj)  —  {aj\bi)  is  the  antisymmetrized  two-electron  integral, 
and  f  -  is  the  general  Fock  matrix  element.  Our  formulation  is  completely  general  [21], 
admitting  non-SCF  orbitals  where  ff,  or  may  be  nonzero.  This  generality  has 
found  a  useful  implementation  in  coupled-cluster  calculations  using  localized  orbitals 
[10]. 

A  consideration  of  the  expansion  of  the  operator  exp(X|  +  T2)  demonstrates  the 
inclusion  of  all  single  and  double  excitations  and  the  disconnected  contributions  of 
triple,  quadruple,  and  higher  excitations.  In  the  quadruple  excitation  case,  the  dis¬ 
connected  contribution  V2  Xt|0o)  is  clearly  the  dominant  part,  making  the  neglect 
of  X4  a  good  approximation  [3,27],  Based  upon  perturbation  theory  arguments,  with 
SCForbitais  the  connected  triple  excitation  parts  that  would  arise  from  Xj  are  expected 
to  be  larger  than  those  from  the  disconnected  X(X2  and  Xf  contributions  [?S], 
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Projection  of  ^ccso  onto  just  single  and  double  excitations  decouples  the  CCSD 
equations  to  terms  quartic  in  T\,  quadratic  in  T2,  and  the  cubic  term  T\Tj-  However, 
the  wave  function  is  not  truncated,  retaining  its  exponential  form.  These  nonlinear 
equations  are  solved  iteratively  using  Pade  approximants  [29]  and  a  recently  proposed 
reduced  linear  equations  method  [30]  to  assist  the  convergence.  The  lowest  order  it¬ 
erations  correspond  to  second-  and  third-order  perturbation  theory  and  the  fourth- 
order  model  limited  to  single,  double,  and  quadruple  excitations  SDQ-MBPT(4) 

13]. 


3.  Results  and  Discussion 

The  insertion  of  Be  into  H2  can  be  a  complicated  process.  The  Be  atom  is  well  known 
for  the  quasidegeneracy  of  the  2s  and  2 p  orbitals,  so  it  is  usually  thought  to  have  two 
important  configurations:  the  ls22s2  and  the  ls22p2  set.  The  comparative  weights 
of  these  two  configurations  are  very  much  a  consequence  of  the  orbital  set  used,  with 
natural  orbitals  [31],  e.g.,  weighting  the  second  configuration  higher  than  is  the  case 
with  SCF  orbitals.  However,  only  the p1  components  allow  the  BeH2  bond  to  be  formed, 
according  to  orbital  symmetry  rules. 

The  insertion  of  Be  into  H2  also  requires  that  the  Hi  bond  be  effectively  broken. 
Yet,  at  large  internuclear  separation,  the  Iff*  configuration  for  H2  becomes  nearly 
degenerate  with  the  1  a \  configuration,  causing  both  configurations  to  be  almost  equally 


Figure  1.  Comparison  of  fv-mcscf  and  FCl  for  the  two  lowest  'Wi  states  of  BeH.  in 
Ch,  symmetry,  (v)  FuK  valence  MCSCFl  '/ti;  (□)  full  valence  MCSCF2  ',«i;  (O.  A) 

FCl. 
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Table  I.  Contracted  Gaussian  basis  used  for 
the  ten  orbital  BeHj  model  problem. 


Se  Is 


Is¬ 

is' 

2p 

H  Is 

Is' 


Exponent 

1267.07 

190.356 

43.2959 

12.1442 

3.80923 

1.26847 

5.693380 

1.565630 

0.171855 


5.693330 

1.555630 

0.171353 

19.2406 

2.8992 

0.6534 

0.17760 


Contraction 

Coefficient 

.001940 

.014786 

.071795 

.236348 

.471763 

.355183 

-0.028876 

-0.177565 

1.071630 


.004836 

.144045 

.949692 

.032828 

.231208 

.B17238 

1.000 


0.057181  1.000 


weighted  in  the  ci  description  of  the  molecule.  Consequently,  a  rational  first  ap¬ 
proximation  to  the  computation  of  the  interaction  energy  along  a  reaction  path  for 
BeH^  would  seem  to  require  at  least  the  configurations  obtained  by  coupling  together 
these  four  configurations  of  the  two  components.  In  BeH;,  though,  there  are  even  more 
complications. 

Since  the  perpendicular  insertion  of  Be  into  H;  requires  /7-orbital  participation  on 
Be,  there  is  a  promotion  from  Be(2s2)  to  Be(2p:)  near  the  critical  geometry.  This 
causes  the  principal  configuration  in  BeHx  to  change  from  la*2a*3aj  to 
In  a  single  reference  model,  one  of  these  two  very  important  configurations  must  be 
treated  in  the  complementary  space,  placing  a  severe  burden  on  any  single  reference 
model  to  describe  the  insertion  reaction.  To  emphasize  the  multiconfigurational  aspect 
of  this  problem,  at  the  peak  on  Figure  1 ,  which  corresponds  to  a  Be  to  Hj  distance  of 
2.75  bohr,  configuration  b  (with  SCF  orbitals  from  b)  has  a  coefficient  of  0.724,  while 
the  a  coefficient  is  0.560.  In  addition,  the  configuration  la,  52a,  -3a,  5a also  has  a 
large  coefficient  of  0.239.  When  Be  is  moved  to  3.0  bohr,  a  has  the  larger  coelficieni 
of  0.323  compared  to  0.294  for  b.  A  final  complication  is  that  both  of  these  RHF  ref¬ 
erence  functions  are  unstable  since  there  arc  lower-lying  UHF  solutions. 

The  contracted  Gaussian  basis  set  is  defined  in  Table  I.  The  sampling  path  for  the 
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Table  II.  Coordinates  of  points  along  the  sampling 
path*  for  the  perpendicular  insertion  of  Be  into  H;. 


Pplnt 

Coordinates 

for  Hj  («.  Y.  2)' 

A 

(0.  1 

2.54,  0.0) 

1 

(0,  » 

2.08.  1.0) 

C 

(0.  * 

1.62,  2.0) 

0 

(0.  2 

1.39.  2.5) 

c 

(0.  i 

1.275,  2.75) 

F 

(0.  2 

1.16.  3.0) 

s 

(0.  2 

0.93.  3.5) 

N 

(0.  2 

0.70.  4.0) 

I 

(0.  -• 

0.70.  6.0) 

J 

(0.  2 

0.70,  20.0) 

•  This  sampling  path  is  selected  to  allow  simple  com¬ 
parisons  of  calculated  energies  along  a  path  which  samples 
points  near  the  reaction  path. 
k  Be  is  located  at  (0.0.0.Q.0.0). 


Civ  insertion  energy  surface  is  defined  in  Table  II.  The  r  coordinate  is  the  distance 
of  Be  from  the  midpoint  of  the  H2  bond.  The  y  coordinate  defines  the  location  of  the 
two  H  atoms.  This  sampling  path  has  determined  as  follows.  MCSCF  calculations  were 
performed  to  locate  the  vicinity  of  the  transition  state  for  the  perpendicular  insertion 
of  Be  into  Hi.  On  a  graph  where  the  Be — H2  distance  is  the  z  axis  and  the  H — H 
distance  is  the  y  axis,  a  straight  line  was  drawn  between  the  BeHi  equilibrium  point 
and  the  transition  state  geometry.  A  second  straight  line  was  drawn  parallel  to  the 
z  axis  at  the  Hi  equilibrium  distance  and  the  intersection  with  the  other  line  defines 
the  corner  of  the  sampling  path.  Clearly,  this  is  not  the  reaction  path,  but  has  the 
advantage  of  permitting  a  comparison  of  energies  determined  with  different  models 
at  the  same  point  for  all  methods  while  passing  through  regions  of  the  potential  energy 
surface  near  the  reaction  path.  The  true  reaction  path  for  the  unconstrained  insertion 
of  Be  into  H2  need  not  have  Civ  symmetry. 

Before  discussing  the  CCSD  results,  we  show,  in  Figure  I ,  a  comparison  of  the  full 
CI  with  a  full  valence  (FV-MCSCF)  (also  called  a  complete  active  space  [32]  or  full 
reaction  space  [33])  calculation  for  the  lowest  two  1  ^-4 1  states  of  BeH:  [25,26].  The 
active  space  for  the  FV-MCSCF  includes  the  2a (,  3a  1,  4a  1,  5a  1, 162.  2f>2.  and  16)  or¬ 
bitals,  leading  to  70  configuration  state  functions.  The  wave  function  is  then  optimized 
in  terms  of  the  orbital  coefficients  and  configuration  state  function  mixing  coefficients. 
This  would  appear  to  be  a  systematic  well-defined  approximation  that  exploits  the 
multireference  approach  to  the  fullest,  short  of  adding  the  remaining  ci  functions. 

The  agreement  between  the  FV-MCSCF  and  full  CI  (FCI)  is  quite  good,  although 
there  is  an  energy  difference  varying  from  2.4  to  5.4  mhartree  that  is  noticeable  even 
on  the  50  mhartree  scale  of  Figure  1 .  There  are  also  different  errors  for  different  parts 
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of  the  sampling  path,  with  the  greatest  errors,  surprisingly,  being  near  the  end  points 
for  the  ground  state. 

Alternatively,  we  can  take  a  single,  albeit  poor,  reference  and  see  whether  CCSD 
can  introduce  the  requisite  corrections.  To  do  so,  we  first  have  to  define  the  reference 
function,  and  this  will  itself  introduce  ambiguities. 

Since  this  is  a  singlet  molecule,  we  will  choose  to  use  a  spin-restricted  Hartree-Fock 
determinant  as  the  reference  determinant.  We  thus  avoid  some  of  the  well-known 
problems  associated  with  the  use  of  unrestricted  Hartree-Fock  determinants  to  de¬ 
scribe  bond  formation  and  breaking.  However,  the  use  of  space  and  spin  restricted 
Hartree-Fock  reference  determinants  here  introduces  its  own  symmetry  problems. 
While  the  symmetry  problems  associated  with  UHF  functions  are  well  known,  the  use 
of  restricted  Hartree-Fock  in  the  present  Ciu  insertion  reaction  forces  too  much 
symmetry.  In  particular,  the  RHF  method  allows  the  crossing  of 1  ^4 1  determinants 
laT2ai3<J:  (hereafter  called  a)  and  laf2nfl^5  (called  b)  configurations,  resulting  in 
two  ways  to  obtain  an  energy  for  the  '.-It  state.  This  crossing  is  allowed  at  the  RHF 
level  because  of  the  extra  symmetry  element  imposed  by  separating  electrons  into 
orbitals.  Since  the  full  Cl  is  invariant  to  the  choice  of  molecular  orbitals  and  allows 
both  configurations  to  detet.r'ne  their  weight  in  the  ci  solution,  the  full  Ci  naturally 


Table  III.  Comparison  of  CCSD.  MBIT.  and  FV-MCSCF  with  full  Cl  results  for  the 
perpendicular  insertion  of  Be  into  H,  (C^,  symmetry).* 


Nlil 

sew 

jao-mrri*) 

«* 

n-*c% cr 

U* 

CC» 

o 

TCI 

-is.:*i7 

-IS. 7774 

x.o 

-is.r** 

».o 

•IS. 77*4 

4.1* 

-11.7742 

-IS. 4944 

-IS.7JS* 

1.4 

-IS. 7127 

2.4 

•IS. 73*4 

0.2* 

-IS. * '72 

-1S.424* 

•IS. 4714 

1.0 

-IS. 470* 

2.S 

-IS. *7*) 

4.32 

-IS. 41*4 

-is  .Hit 

•IS.  41 S) 

4.7 

-U.42  93 

1.2 

•  !9.*aJ  7 

0.13 

-IJ.4774 

-is.  mi 

-IS. 540* 

14.1 

-17.7M3 

l.t 

-1S.S97I 

1.4 

*13  ■  Ml* 

-IS. sits 

-IS. SMI 

u.t 

-IS.  SMI 

1.2 

-IS. *411 

-.25 

-1*40  * 

-IS. AMS 

-IS. *714 

14.  S 

-IS.  SMI 

1.2 

-IS. 5  Ml 

1.) 

-IS. *429 

-IS.  SMS 

-11.4274 

7.4 

•IS. 4179 

1.9 

-19- #7-1 

4.  *4 

•11.42S4 

-is.tisr 

-IS. 447! 

J.  7 

-IS.*4M 

2.1 

-IS.**!* 

4.21 

-13.4*12 

-IS. MM 

-IS. 7311 

1.4 

-is.rni 

J.* 

-IS. 7341 

4.12 

-IS.  7  347 

-IS. MS* 

-1S.7S4# 

1.4 

-IS. MM 

1.3 

-IS  **04 

4.01 

-19.7*0* 

-13.4974 

-IS.7S44 

l.J 

-11.7974 

9.7 

-19.741* 

4.41 

-11.7474 

'There  were  1374  configuration  functions  used  in  the  FCI  expansion.  and  141  amplitudes 
in  the  lai2a;?aj  CCSD  expansion.  All  numbers  are  in  a.u. 

*  Energy  difference  from  the  fci  i»  in  kcal/ mol. 

•  MBPT/CCSD  energies  for  points  A  —  D  are  determined  from  the  SCF  reference  configuration 

Idi2djlbi. 

4  Energy  obtained  from  pure  symmetry  SCF  reference  lu;2ailf>; 

•  Energy  obtained  from  a  broken  symmetry  SCF  relerence  configuration.  In  this  case  there 
arc  381  amplitudes. 

1  Energy  obtained  from  pure  svmmetry  SCF  reference  ln;2u;3<i; 

*  MBPT/CCSD  energies  lor  points  F  —  J  are  determined  from  the  SCF  reference  configuration 
Idj2d;3af. 
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Figure  2.  Comparison  of  (a)  SCF  and  (X)  CCSD  results  with  (O)  FCI  for  the  ground 
'A  I  state  of  BeH2  in  Cu  symmetry.  Different  RHF  reference  functions  are  used  at  dif¬ 
ferent  points  on  the  sampling  path.  See  Table  Ill  for  detailed  specifications.  {□)  SCF 
symmetry  broken;  (•)  CCSD  symmetry  broken. 

resolves  this  problem.  MCSCF  does  too,  since  one  set  of  optimum  orbitals  is  obtained 
with  participation  from  several  important  configurations. 

In  any  nonoptimized  orbital  limited  Cl,  though,  this  problem  is  not  resolved.  Unlike 
full  Cl,  a  truncated  ci  is  not  invariant  to  this  kind  of  orbital  modification,  leading  to 
different  solutions  depending  upon  which  reference  function  is  chosen  to  define  the 
orbitals.  Since  CCSD  is  also  a  function  of  its  reference  configuration,  in  spite  of  many 
other  attractive  invariance  properties  [21  ],  it  is  also  possible  to  obtain  separate  CCSD 
results  for  each  choice  of  reference  configuration. 

In  Table  III  and  Figure  2,  we  use  the  reference  function  b  for  points  A-D  and  a 
for  points  F-J ,  whose  geometries  are  defined  in  Table  I.  Besides  point  E,  which  will 
be  discussed  separately,  this  always  corresponds  to  the  lowest  CCSD  energy.  This  also 
reflects  the  lowest  RHF  energy  except  at  point  D ,  where  the  a  configuration  SCF  energy 
is  lower  by  0.126  hartree  while  the  CCSD  result  with  the  b  configuration  is  slightly 
lower  by  0.2  mhartree. 

As  illustrated  in  Figure  2  and  Table  III,  the  CCSD  results  are  very  close  to  the  full 
Cl,  differing  by  less  than  1  keal/mol  except  at  point  E.  This  is  a  surprisingly  good  result 
for  a  single  reference  approach  in  a  problem  with  such  obvious  multiconfigurational 
character.  CCSD  has  141  coefficients  compared  to  1574  configurations  in  the  full  Cl. 
Furthermore,  the  CCSD  energies  are  significantly  better  than  the  FV-MCSCF  energies, 
the  latter  of  which  results  from  taking  full  advantage  of  the  multiconfigurational 
character  within  its  70  configuration  reference  space.  Previous  calculations  [25]  in- 
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Abstract 

A  series  of  Many-Body  Perturbation  Theory  (MBPT)  and  Coupled-Cluster 
(CC)  calculations  were  performed  on  the  ethylene  molecule  employing  both 
the  standard  canonical  SCF  and  localized  orbitals  (LO) .  The  geometry  and 
(4s2pld/2slp)  basis  set  are  Identical  to  that  used  in  the  over  one  million  con¬ 
figuration  "bench  mark"  Cl  calculation  of  Saxe  et.  al.  Using  SCF  orbitals 
we  computed  full  4th  order  MBPT,  (i.e.  SDTQ-MBPT  (4)),  CC  doubles  (CCD)  and 
CC  singles  and  doubles  (CCSD)  energies  and  compare  with  their  Cl  results. 

In  addition,  LOs  were  generated  in  a  manner  similar  to  the  semi-empirical 
PCILO  procedure  and  SDQ-MBPT(4)  and  CCSD  calculations  were  carried  out 
within  this  basis.  Though  the  SCF  and  LO  reference  determinant  energies 
differ  by  0.29706  hartrees,  the  CCSD  energy  difference  is  just  0.00170 
hartrees.  Comparing  with  Saxe  et.  al.  our  most  extensive  SCF  orbital  calcu¬ 
lation,  CCSD  plus  4th  order  triples,  was  found  to  be  lower  in  energy  than 
the  Cl  result  by  -0.00533  hartrees. 
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Introduction 


Orbitals  localized  on  one  or  two  atomic  centers  possess  a  number  of 
advantages  over  the  typical  delocalized  canonical  Self-Cons istent-Fleld 
(SCF)  orbitals  for  the  study  of  molecular  potential  energy  hypersur- 
faces[l].  They  are  conceptual;  they  are  concentrated  largely  on  a  single 
atom  or  two  nearest  neighbor  atoms;  and,  they  approximate  traditional  che¬ 
mical  bonding  concepts.  For  example,  within  an  LO  basis  we  might  hope  to 
develop  bond  additivity  schemes  that  would  apply  for  molecular  properties 
such  as  the  total  energy  and  the  dipole  moment.  Transferability  of  calcu¬ 
lated  properties  of  molecular  fragments  between  different  molecules  might 
also  be  expected  to  be  more  accurate  In  a  localized  basis.  Lastly,  but 
most  important  from  a  computational  polnt-of-view,  is  the  reduction  In 
the  number  of  molecular  orbital  (MO)  integrals  required  in  a  correlated 
energy  calculation.  In  a  typical  delocalized  SCF  basis  the  number  of  non¬ 
zero  integrals  grows  as  the  fifth  power  of  n  where  n  Is  the  number  of  MOs 
Included.  However,  in  a  localized  basis  the  number  of  MO  integrals 
increases  only  linearly  with  the  number  of  bonds  for  extended  systems[2]. 
This  reduces  the  required  time  for  both  the  Integral  transformation  and  the 
MBPT/CC  calculation  substantially. 

LOs  can  be  generated  In  a  variety  of  ways,  such  as  the  energy  Invariant 
transformations  among  SCF  orbitals  as  in  the  procedures  of  Foster  and  Boys 
[3]  and  of  Edminston  and  Ruedenberg  [4].  From  the  viewpoint  of  pertur¬ 
bation  theory,  these  procedures  have  the  advantage  that  the  single 
reference  energy  Is  optimum,  suggesting  better  approximations  In  low-order 
perturbation  theory.  Also,  the  prescription  for  forming  the  locallzed- 
orbitals  is  independent  of  geometrical  considerations,  providing  a  well- 
defined  procedure  when  following  a  bond  formation,  for  example.  However, 


the  SCF  procedure  is  more  complex  and  requires  the  generation  of  additional 
molecular  integrals.  From  an  infinite  order  viewpoint,  models  like  con¬ 
figuration  interaction  singles  and  doubles  (CISD)  and  coupled  cluster 
singles  and  doubles  (CCSD)  are  invariant  to  these  localized  orbital  trans¬ 
formations,  since  only  occupied  or  excited  orbitals  are  mixed  among  them¬ 
selves,  so  the  final  energies  will  be  the  same  as  when  using  SCF  orbitals 
provided  the  MO  basis  is  not  truncated. 

A  more  primitive  category  of  localized  orbitals  are  LOs  constructed  to 
directly  represent  the  chemical  bonds  in  a  molecule  as  in  the  PCILO  proce¬ 
dure!  5  J.  These  methods,  built  on  hybrid  and  bond  orbitals,  are  usually 
computationally  quite  simple.  No  SCF  calculation  is  required,  however, 
unlike  localized  SCF  orbitals,  the  CCSD  or  CISD  energy  will  not  be  the  same 
as  in  the  SCF  orbital  case.  In  this  approach  a  minimum  basis  atomic  orbi¬ 
tal  set  is  used  to  form  the  LOs,  hence,  even  if  a  large  basis  (DZP  or 
better)  is  used  to  define  the  AOs,  the  resultant  single  determinant  energy 
does  not  benefit  at  all  from  the  polarization  functions  in  the  basis.  It 
has  been  shown  that  this  energy  difference  is  of  the  same  order  of  magni¬ 
tude  as  the  correlation  energy  Itself [2],  thereby  placing  a  greater  burden 
on  the  procedure  for  incorporating  corrections  to  the  single-determinant 
reference  function.  However,  we  found  In  a  preliminary  study  using  these 
simple  LO’s  for  CH^  and  ^0  that  If  extensive  correlated  calculations  are 
performed,  such  as  CCSD,  very  good  agreement  with  analagous  calculations 
employing  SCF  orbitals  is  obtained  even  including  an  accurate  description  of 
bond  stretching[2] .  In  the  current  study  we  extend  this  treatment  to  CjH^. 

Correlation  effects  in  ethylene  beyond  the  CISD  level  have  recently 
been  examined  by  Saxe,  Fox,  Schaefer  and  Handy  (SFSH)[6J,  using  a 
(4s2pld/2slp)  contracted  Gaussian  basis  set.  In  addition  to  all  single  and 
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double  excitations  from  the  *A  reference,  a  substantial  subset  of  all 

S 

triple  and  quadruple  excitations  were  Included  In  a  massive  1,046,758  con¬ 
figuration  Cl  calculation.  This  set  should  include  the  most  Important 
triple  and  quadruple  excitations  and  should  give  a  very  good  estimate  to 
both  the  Cl  containing  1  though  4  fold  excitations  (CISDTQ)  and  the  full  Cl 
results.  Since  this  Cl  provides  some  values  for  higher  categories  of  exci¬ 
tation,  some  interesting  comparisons  with  MBPT/CC  results  are  possible, 
which  we  report. 

Comparison  of  MBPT/CC  Results  With  Cl 

In  our  study  of  ethylene  we  employed  the  same  basis  set  and  geometry  as 
SFSH[6].  The  contracted  Guasslan  set  used  is  the  standard  Huzlnaga-Dunning[7] 
double  zeta  basis  set  augmented  with  a  set  of  d  functions  on  each  carbon 
and  a  set  of  p  functions  on  each  hydrogen.  Precisely,  the  DZP  basis  sets 
for  carbon  and  hydrogen  are  designated  (9s5pld/4s2pld)  and  (4slp/2slp) 
respectively.  A  factor  of  1.2  was  used  to  scale  the  DZ  hydrogen  basis 
functions.  Six  cartesian  d  functions  were  placed  at  each  carbon  center 
and  each  has  an  exponent  of  0.75.  The  hydrogen  p  exponent  is  1.0.  The 
geometry  of  ethylene  was  chosen  to  be  r(C-C)  =  1.330A,  r(C-H)  ■  1.076A, 

<HCH  “  116.6  which  reproduces  SFSH's  nuclear  repulsion  energy  of 
33.51358956  to  all  reported  digits[6]. 

Using  this  basis  set  and  geometry  a  number  of  energy  calculations  using 
LO’s  and  SC F  orbitals  were  performed  on  the  1A^  ground  state  of  All 


calculations  involving  the  SCF  orbitals  were  done  In  full  D2h  symmetry  and 
correspond  to  the  configuration. 


The  SCF  energy  is  computed  to  be  -78.050530  hartrees  which  agrees  with  the 
reported  value  of  SFSH.  Beyond  the  SCF  level  we  computed  the  full  4th 
order  MBPT  (designated  SDTQ-MBPT(4)  where  SDTQ  represents  all  single, 
double,  triple  and  quadruple  excitation  diagram  contributions) [8,9],  coupled- 
cluster  doubles  (CCD) [8]  and  CCSD[10]  energies  using  the  SCF  orbitals.  (See 
Reference  [11]  for  the  procedure  for  evaluating  the  triple  excitation 
diagrams).  Consistent  with  SFSH,  all  of  our  correlated  calculations  kept 
the  two  carbon  Is  core  orbitals  “frozen".  These  results  are  presented  in 
Table  1,  along  with  selected  Cl  results  of  SFSH. 

A  short  summary  about  the  SFSH  calculation  is  in  order.  Initially 
they  carried  out  a  standard  one  reference  Cl  including  all  single  and 
double  excitations  (CISD)  in  the  SCF  orbital  basis.  One  natural  orbital 
(NO)  iteration  was  next  performed  in  this  5057  configuration  space  to 
obtain  a  "better"  set  of  orbitals  for  the  remaining  Cl  calculations.  To  go 
beyond  CISD  the  orbital  space  was  partitioned  into  two  groups:  the 
valence  space  and  the  remainder.  Their  valence  space  is  composed  of  the 
occupied  orbitals  plus  the  l^g*  ^blu*  ^b2u*  ^ag’  ^3g  and  ^blu  virtual 
orbitals  and  corresponds  to  the  minimum  basis  subspace  with  the  core  fro¬ 
zen.  This  partitioning  is  unsatisfactory  for  SCF  orbitals  since  the  gap 
between  the  orbital  energies  of  the  highest  valence  and  lowest  non-valence 
orbital  is  quite  small.  This  implies  that  certain  energetically  close 
orbitals  will  be  treated  in  an  unequal  manner.  Using  the  natural  orbitals 
Instead,  yields  a  larger  separation  and  a  more  satisfactory  partitioning. 

Two  very  large  Cl  calculations  were  carried  out  using  this  par¬ 
titioning.  In  the  smaller  calculation  the  primary  natural  configuration 
and  all  valence  space  single  excitations  were  used  as  reference  functions 
(this  reference  calculation  is  termed  CIS)  and  then  all  singles  and  doubles 
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from  the  reference  functions  were  used  in  a  Cl  (CISfSD).  This  procedure 
yields  109,473  configurations  including  all  single,  double  and  a  subset  of 
triple  exciations  with  respect  to  the  NO  reference.  In  their  largest  Cl 
all  single  and  double  excitations  in  their  valence  space  plus  the  primary 
configuration  were  selected  as  reference  functions  and  all  single  and 
double  excitations  were  generated  out  of  this  subspace  (CISDfSD).  This 
1,046,768  configuration  basis  is  conposed  of  all  single  and  double 
excitations  and  a  subset  of  all  triple  and  quadruple  exciations  from  the  NO 
reference. 

In  their  preliminary  CISD  calculations  a  correlation  energy  of  -0.27758 
was  found  using  SCF  orbitals,  and  -0.27748  with  NO  orbitals.  Comparing  our 
MBPT/CC  results,  these  values  are  close  to  our  second  order  MBPT  energy, 

E2,  of  -0.27222  hartrees.  Adding  in  Eg,  the  third  order  contribution,  a 
total  correlation  energy  of  -0.29655  hartrees  i3  obtained.  In  their  most 
complete  calculation,  CISDPSD,  their  total  energy  is  -78.35451  hartrees 
yielding  a  correlation  energy  of  0.30398  hartrees.  From  our  fourth-order 
calculation,  SDQ-MBPT(4)  and  infinite  order  CCSD  result,  we  compute  corre¬ 
lation  energies  of  0.30004  and  0.30116  hartrees,  respectively,  or  about 
98.7Z  and  99.lt  of  the  CISDK3D  correlation  energy. 

In  CCSD  and  Its  fourth-order  approximation,  SDQ-MBPT(4)  [8],  the 
contribution  from  the  triple  excitation  diagrams  is  neglected.  With  SCF 
orbitals  these  diagrams  first  appear  in  MBPT  in  the  4th  order  (see  Ref. 
[8,9])  and  we  designate  this  energy  contribution,  E^.  Computing  E^  we 
obtain  an  energy  of  -0.00815  hartrees.  This  is  effectively  identical 
with  the  E^t  contribution  of  Frisch  et.  al.  of  -0.00815  hartrees  in  the 
slightly  smaller  6-31G*  basis  set  [12].  Substracting  the  CIS+SD  and  CISD 
energies  yields  an  estimate  of  the  importance  of  triple  excitations  in  the 


Cl  framework  of  -0.0070L  hartrees,  which  Is  similar  to  our  MBPT  result. 

Adding  E^t  to  the  SDQ-MBPT(4)  and  CCSD  energies  our  two  best  estimates  of 
the  correlation  energy  are  -0.30840  and  -0.30931  hartrees  respectively.  The 
corresponding  total  energies  for  the  SDTQ-MBPT(4)  and  CCSD+T(4)  models  are 
-78.35873  and  -78.35984  hartrees.  These  energies  are  below  the  CISD+SD  value 
by  4.4  and  5.3  mh,  respectively.  Since  the  triple  excitation  selection  can 
account  for  about  1  mh  of  this  difference,  it  would  appear  that  as  much  as 
3-4  mh  arises  from  the  subset  of  quadruple  excitation  configurations  chosen. 

To  analyze  this  in  a  bit  more  detail,  we  can  make  some  approximate  com¬ 
parisons  for  the  contribution  of  quadruple  excitations  to  the  Cl.  As  we 
have  discussed  in  detail  elsewhere[13] ,  the  contribution  of  quadruple  exci¬ 
tations  to  Cl  is  approximately  given  by  the  difference  of  the  fourth-order 
linked  and  unlinked  quadruple  excitation  diagrams  of 

Eci(Q)«E4q«2A  (1) 

The  value  of  the  unlinked  daigrams  depends  upon  A®<<|>^  |  <|>^>,  which  is  the 
overlap  of  the  first-order  perturbed  wavefunction.  This  formula  forms 
the  basis  for  the  widely  used  Davidson's  approximation  [14]  for  E^j(Q) 

(where  EgjCQJsd^A)  as  we  proved  elsewhere  [15,9].  In  our  calculation, 
E2A*-0.02647  hartrees.  Using  E^q  from  Table  1,  we  obtain  ECI(Q)»-20.8  mh.  An 
appropriate  estimate  for  EC^(Q)  from  the  Cl  calculation  is  given  by  the  dif¬ 
ference  between  CISD+5D(N0)  and  CIS+SD(N0)  which  is  -19.5  mh,  which  is  in 
good  agreement.  Although  this  number  arises  largely  due  to  the  addition  of 
quadruple  excitations  into  the  Cl,  some  additional  triple  excitation  terms 
are  also  included.  Again  using  an  estimate  ~1  mh  for  the  additional  triple 
excitations  in  CISD+SD  gives  ~2-3mh  as  the  energy  lost  in  the  quadruple 
selection,  which  is  consistent  with  our  above  observation.  Clearly,  any 
such  estimates  mixing  infinite-order  and  finite-order  calculations  as  well 


as  different  (but  similar  SCF  and  natural)  orbitals  are  not  rigorous,  but 
appear  to  be  suggestive  of  the  size  of  the  selection  error. 

The  over  one  million  configuration  CISD+SD  calculation  of  SFSH  was  per¬ 
formed  primarily  as  a  demonstration  of  the  feasibility  of  such  large  com¬ 
putations  and  used  code  largely  designed  for  just  such  a  purpose,  hence 
____Ximi«g  comparisons  are  not  entirely  appropriate.  Recognizing  this,  we  can 
make  a  few  very  approximate  comparisons.  The  CISD+-SD  calculation  was 
carried  out  on  a  Harris  800  miniconputer  and  required  just  over  100  hours 
of  CPU  time  (13  hours  per  iteration).  The  SDQ-MBPT(4)  and  CCSD  timings  by 
conparlson  took  0.2  and  2.2  hours,  respectively,  on  an  IBM  3033N  computer. 
The  IBM  runs  electronic  structure  calculations  at  between  2  to  5  times  the 
speed  of  the  Harris  model  thus  suggesting  about  20  hours  for  the  Cl  on  the 
3033N,  or  roughly  100  times  as  long  as  the  SDQ-MBPT(4)  and  about  10  times 
as  long  as  the  CCSD  calculations.  Computing  the  E^  energy  is  inherently 
harder  than  In  the  SDQ-MBPT(4)  or  CCSD  cases  since  finding  E^T  is  propor¬ 
tional  to  n^  while  the  latter  are  proportional  to  n^  operations  where  n 
is  the  number  of  orbitals.  Even  so,  E^T  required  almost  1.0  hour  to  com¬ 
pute.  Our  best  calculations,  SDTQ-MBPT(4)  and  CCSD+T(4)  therefore,  reouire 
1.2  and  3.2  hours,  respectively,  or  6.0/E  and  16. OX  of  the  CISD+SD  time. 

The  time  savings  occurs  essentially  due  to  the  fact  that  only  a  number  of 
coefficients  equal  to  all  single  and  double  excitations  are  required  in 
CCSD[10],  although  non-linear  terms  are  considered,  compared  to  the  over  1 
million  coefficients  in  the  Cl.  CISD  has  a  similarly  small  number  of  coef¬ 
ficients  and  requires  less  computer  time  than  the  MBPT/CC  calculations. 


—  8  — 


Comparison  of  Localized  Orbital  and  SCF  Orbital  Results 

Our  second  series  of  calculations  on  ethylene  used  simple  bond- 
antlbond  localized  orbitals.  These  LOs  were  constructed  as  follows.  (A 
more  detailed  discussion  is  in  ref.  [2].)  First,  the  2s,  2px  and  2pz  ato¬ 
mic  functions  on  each  carbon  atom  were  hybridized  to  form  a  set  of 
2 

sp  functions  directed  at  the  opposing  carbon  and  the  two  nearest  hydrogen 
atoms.  Next,  linear  combinations  of  hybrid  orbitals  were  formed  to  yield 
bonding  and  antibonding  orbitals.  Precisely,  the  non-orthogonal  C-H,  sigma 
(C-C)  and  pi(C-C)  bonds  and  antibonds  are  [C(sp2 )±H(ls)  ], 

Cc1(sp2)±C2(sp2)]  and  tcl(Py^±C2^Y^ respectively,  which  forces 
equal  weighing  even  for  the  CH  bond.  Non-polar  bonds  are  not  necessary  but 
chosen  purely  for  convenience  in  the  present  calculation.  (See  Ref.  [16] 
for  the  effect  of  non-equal  weights  in  the  PCILO  scheme  and  Ref.  [17]  for 
an  ab-lnltio  study.)  All  polarization  functions  in  the  DZP  basis  are 
Introduced  into  the  calculation  as  excited  orbitals  in  addition  to  the 
anti-bonding  orbitals. 

In  the  PCILO  method  [5],  the  zero-differential  overlap  approximation  is 
normally  used,  however  an  ab-lnltlo,  non-mininum  basis  calculation  requires 
an  orthogonalization  procedure.  Instead  of  carrying  out  a  Lowdin  symmetric 
orthonormalization  [18]  on  the  complete  localized  set,  the  basis  is  first 
divided  into  three  subsets.  These  are  the  core  plus  bonding  orbitals,  the 
antibonding  orbitals,  and  the  remaining  excited  orbitals.  The  third  cate¬ 
gory  is  necessary  to  accommodate  larger  than  minimum  AO  basis  sets.  (The  term 
"excited"  orbitals  refers  to  antibonding  orbitals  plus  virtual  orbitals.) 

Next  each  subset  is  individually  symmetrically  orthonormalized  and  the  anti- 
bonding  block  is  Schmidt  orthonormalized  to  the  bonding  orbitals.  Lastly,  the 
virtual  orbitals  are  Schmidt  orthonormalized  to  both  the  bonding  and  anti- 
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bonding  subsets.  There  are  several  advantages  for  employing  this  orthonor- 
malizlng  procedure.  Perhaps  the  most  Important  is  that  the  bonding  orbitals 
contain  no  contamination  from  the  excited  orbitals.  This  feature  should  faci¬ 
litate  potential  transferability  among  different  calculations.  In  addition, 
since  the  CCSD  energy  is  invariant  to  occupied-occupied  and  excited-excited 
orbital  rotations,  use  of  either  localized  or  delocalized  virtual  orbitals 
would  result  in  energetically  equivalent  results. 

The  reference  energy  for  the  single  determinant  formed  by  doubly 
occupying  the  two  Is  core  and  bonding  orbitals  is  -77.75347  hartrees. 

Within  this  LO  basis  a  number  of  correlated  calculations  are  performed, 
including  SDQ-MBPT(4)  and  CCSD.  As  in  the  calculations  employing  SCF  orbi¬ 
tals,  the  two  core  orbitals  are  kept  frozen. 

Comparing  first  the  SCF  and  LO  reference  determinant  energies  a  difference 
of  -0.29706  hartrees  is  found.  At  first  this  difference  seems  large  since  the 
LOs  should  lead  to  a  reasonable  guess  at  the  reference  determinant.  However, 
since  our  LOs  are  formed  from  linear  combinations  of  the  valence  atomic  func¬ 
tions,  the  polarization  functions  in  the  DZP  set  contribute  nothing  to  the 
reference  determinant  energy,  nor  is  there  any  opportunity  for  the  conparative 
weights  in  the  bond  orbitals  to  change  to  more  accurately  reflect  the  correct 
charge  distribution.  Both  of  these  features  are  introduced  into  the  calcula¬ 
tion  via  the  single  excitation  term  in  the  CCSD  wavefunction[ 10]  (i.e. 
a  t 

The  CCSD  wavefunction, 

’|,CCSD*exP(Tl+r2)  I  V  "  exp(T2)exp(T1)  |  *o>  (2) 

benefits  from  the  fact  that  any  single  determinant  wavefunction,  can  be 
interchanged  with  any  other  via  Thouless’  relation  [19], 


I  ♦0>-«xp(T1)  |  $Q> 


(3) 


Hence  Che  difference  In  energy  (AE  ,  between  A.  ,  *A  and  a  localized 

pol  TLocal  To 

orbital  representation  of  the  SCF  solution,  4>gCp,  Is  largely  taken  Into 
account  by 

4Epol-<*o  I  I  + j  <**  II  ab><VbJ~Vj)  <4> 

l<b 

for  the  converged  amplitudes,  {t^},  where  f  is  the  off-dlagonal  Fock 
operatorflO].  In  the  lowest-order  Iteration,  tj*f  hence  we  can  observe 
the  low-order  contributions  to  AE^^  by  conventional  perturbation  theory 
arguments.  The  second-order  contributions  to  AEpQ^  is  -0.23008  hartrees 
giving  a  comparison  value  of  -77.98355  hartrees,  which  is  0.06698  hartrees 
above  ES(,p.  At  convergence,  this  is  inproved  to  -78.04475,  or  0.00578 
hartrees  ibout  ESCp. 

Our  best  estimate  of  the  correlation  energy,  the  SCF  orbital  CCSD+E^-j, 
value  of  -0.30931  hartrees,  is  only  marginally  larger  than  the  SCF-L0 
reference  determinant  energy  separation.  Therefore,  using  this  particular 
LO  set,  the  CCSD  procedure  recovers  twice  as  much  energy  as  in  the  analo¬ 
gous  SCF  case.  However,  the  rates  of  convergence  of  the  LO  and  SCF  orbital 
MBPT  calculations  should  be  different.  Comparing  E2  first,  we  find  that  the  L0 
value  of  -0.48326  hartrees  is  substantially  larger  than  the  SCF  orbital  value 
of  -0.27222  hartrees,  but  represents  81.0?  of  the  CCSD  energy  correction  compared 
to  90.4%  in  the  SCF  orbital  basis.  Similarly,  third-order  LO  improves  this  to 
93.6%  compared  to  98.5%  for  the  SCF  case.  Adding  in  the  fourth-order  terms 
due  to  single,  double,  and  quadruple  excitations  obtains  97.7%  in  the  LO  case 
and  99.6%  for  SCF  orbitals,  still  causing  a  discrepancy  of  14  mh  to  l  mh, 
respectively.  However,  as  seen  in  Table  2,  at  convergence  the  CCSD  energies 
differ  by  only  1.7  mh  regardless  of  orbital  basis  set. 
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In  CjH^  all  bond  orbitals  are  chosen  to  be  equally  weighted  despite  the 
electronegativity  differences  between  C  and  H.  We  have  made  a  study  of  how  the 
CCSD  results  with  LO's  change  as  a  function  of  the  bond  orbital  weighting 
factors,  which  will  be  presented  elsewhere  [17].  The  effect  of  triple  excitations 
with  localized  orbitals  requires  the  computation  of  non-Hartree-Fock  triple 
excitation  diagrams.  We  have  presented  the  theory  and  procedure  for  this 
calculatlon[ll]  and  are  currently  developing  the  computer  code. 

Summary 

In  conclusion,  even  with  very  crude  localized  orbital  descriptions  of 
C2H4,  CCSD  is  able  to  recover  about  the  same  energy  correction  as  when 
using  SCF  orbitals.  This  is  primarily  due  to  the  effect  of  T^  in  the 
coupled  cluster  ansatz,  which  is  able  to  introduce  the  requisite  polariza¬ 
tion  effects  into  the  calculation.  Clearly,  with  Foster-Boys[ 3]  or 
Edmistoir-Ruedenberg[4]  LO's,  much  better  low  order  MBPT  approximations  would  be 
possible. 

From  comparisons  between  MBPT/CC  and  the  >106  configuration  Cl  study  of 
C2H4  by  SFHS,  we  show  that  CCSD  and  full  fourth-order  perturbation  theory, 
SDTQ-MBPT(4) ,  give  quite  good  results.  From  the  SDTQ-MBPT(4)  and  CCSD  value, 
it  is  possible  to  suggest  that  configuration  selection  among  the  triple  and 
quadruple  excitations  as  accomplished  by  SFHS,  results  in  only  a  small  2-5  mh 


error  for  this  molecule. 
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Table  1 

C-H, :  Summary  of  SCF  Orbital  MBPT/CC  Energy  Calculations  Including  the  NO 

Cl  Calculations  of  Saxe,  Fox,  Schaefer  and  Handy  [6].  (DZP  Basis  See) 


Type  of  Calculation 

Correlation  Energy 

Cumulative  Total 

E 

-0.0 

-78.05053 

SCF 

E  (NO) 

-78.04942 

Ref 

CISD(SCF) 

-0.27758 

-78.32810 

CISfSD(NO) 

-0.28448 

-78.33501 

CISD+SD(NO) 

-0.30398 

-78.35451 

E2 

-0.27222 

-78.32275 

e3 

-0.02433 

-78.34708 

E4S 

-0.00163 

E4D 

-0.00751 

E4T 

-0.00815 

E4Q 

0.00563 

SDQ-MBPT(4) 

-0.30004 

-78.35057 

SDTQ-MBPT(4) 

-0.30840 

-78.35873 

CCD 

-0.29919 

-78.34972 

CCSD 

-0.30116 

-78.35169 

ccsd+e4t 

-0.30931 

-78.35984 

Table  2 


C2H4:  Localized  Orbital  MBPT/CCM  Results  (DZP  Basis  Set) 


Type  of  Calculation 

Energy  Correction 

Cumulative  Total  Energy 

E 

-0.0 

-78.05053 

SCF 

E 

-77.75347 

Ref 

E2 

-0.48326 

-78.23673 

E3 

-0.07531 

-78.31204 

E4(SDQ) 

-0.02409 

-73.33612 

SDQ-MBPT(4) 

-0.58266 

-78.33612 

CCSD 

-0.59652 

-78.34999 

